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We calculate the Raman response for the Kitaev spin model on the H-Q, H-l, and H-oo harmonic 
honeycomb lattices. We identify several quantitative features in the Raman spectrum that are 
characteristic of the spin liquid phase. Unlike the dynamical structure factor, which probes both the 
Majorana spinons and flux excitations that emerge from spin fractionalization, the Raman spectrum 
in the Kitaev models directly probes a density of states of pairs of fractional, dispersing Majorana 
spinons. As a consequence, the Raman spectrum in all these models is gapless for sufficiently 
isotropic couplings, with a low-energy power law that results from the Fermi lines (or points) of the 
dispersing Majorana spinons. We show that the polarization dependence of the Raman spectrum 
contains crucial information about the symmetry of the ground state. We also discuss to what 
extent the features of the Raman response that we find reflect generic properties of the spin liquid 
phase, and comment on their possible relevance to a—, ft—, and 7 —LUIrOa compounds. 


I. INTRODUCTION 

Quantum spin liquids (QSLs) have been a topic of 
significant interest for four decades since they were in¬ 
troduced in the seminal work of Anderson.^ The search 
for these exotic phases of matter has focused on low¬ 
dimensional and highly frustrated magnets, where con¬ 
ventional ordering is suppressed and quantum fluctua¬ 
tions are large, opening up the possibility of a ground 
state that preserves all the symmetries of the underly¬ 
ing spin Hamiltonian. In lieu of a conventional magnetic 
order associated with broken symmetry, these symmetry¬ 
preserving QSL phases are characterised by a combina¬ 
tion of properties known as topological order?^'^ 

Recent years have seen remarkable progress both in 
the theoretical understanding of the nature of QSLs and 
in identifying realistic Heisenberg models which might 
host them.^“^^ In parallel, a long experimental quest has 
identified a number of materials as spin liquid candidates, 
and evidence for spin liquid physics exists in several two- 
and three-dimensional geometrically frustrated magnets, 
including triangular,kagome,^^ hyperkagome^® and 
pyrochlore antiferromagnets.^®’^® 

In spite of this progress, however, an unambiguous 
identification of QSL phases in these systems remains 
elusive. The difficulty in experimentally identifying QSLs 
comes from the fact that, unlike states with broken sym¬ 
metry, the topological order characteristic of QSLs does 
not admit a local order parameter. This makes it chal¬ 
lenging to identify “smoking gun” experimental signa¬ 
tures. 

However, one distinctive feature of topological order 
that can be probed by conventional methods,such 
as inelastic neutron^^’^^^^® or Raman scattering^’^^^^, is 
fractionalization. In quantum spin liquids there are 
elementary (chargeless) excitations carrying fractional 
quantum numbers relative to the local constituent de¬ 
grees of freedom (charged electrons and spin 1 magnons). 


Therefore only multiple quasiparticles can couple to ex¬ 
ternal scattering probes. With access to the full spec¬ 
trum, one can in principle disentangle the contribu¬ 
tions of individual quasiparticles with different ener¬ 
gies from the resulting multi-particle continuum, and 
make a quantitative comparison between theory and 
experiment. This is in contrast with thermodynamic 
measurements,which only probe the asymptotic 
low energy response. 

The power of such a quantitative comparison was 
strikingly demonstrated in the one-dimensional sys¬ 
tems CUSO 4 • 5 D 2 O and KCuFa.®^^®® Suggestively, the 
experiments®^’®^ see broad spectra characteristic of frac¬ 
tionalized excitations. The fractionalization in these ID 
systems was definitively identified by excellent quantita¬ 
tive agreement between experiments®®’®^ and exact calcu¬ 
lations of the two- and four-spinon dynamical structure 
factor based on an exact solution of the Heisenberg model 
in one dimension using Bethe Ansatz,®®’®® which shows 
fractionalized spin excitations.®® 

Broad spectra have also been predicted and ob¬ 
served in two- and three-dimensional candidate QSL 
materials.However, the level of quan¬ 
titative agreement between experiment and theory 
demonstrated in ID has yet to be achieved in higher 
dimensions.®*®’^® In large part, this is because obtaining a 
reliable calculation of the spin response functions for spin 
liquid states in Heisenberg models has proven difficult 
due to a lack of controlled methods that can be used to 
treat these highly frustrated systems. For the Heisenberg 
QSL candidates theoretical calculations are restricted ei¬ 
ther to numerics, for which obtaining entire spectra in 2D 
and 3D models is quite challenging (see Ref. 9 and refer¬ 
ences therein), or to variational-type ansatz such as slave 
particle®^^®® or resonant-valence-bond treatments.®®’^®® 

However, for a particular class of spin-exchange Hamil¬ 
tonians based on the Kitaev honeycomb model,®® the the¬ 
oretical situation is much more tractable. For these (in- 
tegrable) models both the dynamical structure factor®® 
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and the Raman response^^ can be obtained analytically. 
This is highly significant, since any quantitative features 
of the resulting spectrum that are characteristic of the 
spin-liquid phase can provide a reliable basis for compar¬ 
ison with experiments. 

In this context, an exciting current development in the 
search for spin liquids has been the synthesis of transi¬ 
tion metal compounds belonging to the A2lr03 iridate 
family,^^^®^ and more recently, a-RuCla.®®’^® These ma¬ 
terials form trivalent lattice structures in the harmonic 
honeycomb family and have magnetic moments aris¬ 
ing from Ir'^’*' or Ru^^ ions, which due to strong spin- 
orbit coupling are characterized by Jeff = 1/2 states. 
Edge-sharing IrOg or RuClg octahedra provide 90 ° paths 
for a Kitaev-like super-exchange coupling among mag¬ 
netic moments.Although the known candidate com¬ 
pounds A2lr03, as well as a-RuCl3, actually order at low 
temperature,®^ the presence of a large Kitaev term sug¬ 
gests that these ordered ground states are proximate to 
spin liquid phases - a statement which is also supported 
by recent experiments. 

The tractability of the Kitaev model allows for explicit 
computations of the dynamics.In particular, spins 
fractionalize into two types of degrees of freedom in the 
Kitaev model: dispersing Majorana spinons and gapped 
flux loops. In the unperturbed Kitaev model the flux is 
conserved (static), and the ground state can be viewed 
as a band insulator, or metal, of Majorana spinons in the 
flux background which minimizes the total energy. The 
dynamical structure factor exhibits a gap that can be 
understood as the energy to excite a pair of fluxes (a flux 
loop in 3 D) because spin excitations fractionalize into 
both dispersing Majorana spinons and gapped fluxes. 

These predictions for the 2 D model have already been 
tested experimentally and can be understood within the 
fractionalization of the Kitaev model. While the stan¬ 
dard tool for measuring the dynamical structure fac¬ 
tor - neutron scattering - is challenging in iridates be¬ 
cause of strong neutron absorption, neutron scattering 
data have been recently reported for a-RuCl3.^° Both a 
gapped spectrum and a broad continuum were observed 
and are consistent with the theoretical calculation."*^. On 
the other hand, low-energy photons create pairs of Ma¬ 
jorana spinons (no fluxes), allowing the Raman opera¬ 
tor to directly probe a two-Majorana spinon density of 
states (DOS). Indeed, Raman scattering experiments in 
Q;-RuCl3®^ show a broad scattering continuum consistent 
with this prediction, as well as the prediction of no de¬ 
pendence on polarization.®^’®* Similar but somewhat con¬ 
troversial Raman scattering data have been reported for 
Na2lr03.®® 

The experiments on two dimensional compounds illus¬ 
trate the power of dynamical scattering experiments for 
identifying these elusive phases of matter. Our work fo¬ 
cuses on calculations of dynamical Raman spectra for the 
2 D and 3 D harmonic honeycomb systems, in particular 
on the JL-oo, "H-O, and "H-I lattices realized in the a—, 
/ 3 — and 7—Li2lr03 compounds.®*"’®® We find that the spin 


liquid’s Raman spectrum displays a broad continuum re¬ 
lated to the two-Majorana spinon DOS ( 2 -DOS). In the 
gapless phase there is spectral weight down to zero fre¬ 
quency. In fact, the low energy asymptotic spectrum 
reflects the Majorana spinon Fermi surface topology. In 
addition, the 3 D lattice structures lead to a richer po¬ 
larization dependence as well as more spectral features 
than their 2 D counterpart, due to the larger unit cells 
and lower symmetry. We show that most of these fea¬ 
tures result from the high symmetry of the spin liquid 
phase, together with the particular form of the Raman 
operator in systems with only bilinear spin-exchange in¬ 
teractions, and are therefore expected to be characteristic 
of the QSL phase in these systems. 

The structure of the paper is as follows. In section II, 
we review the Kitaev model on the harmonic honeycomb 
lattices and give some relevant details of their structures. 
Section III reviews the theory of Raman scattering in Ki¬ 
taev spin liquids. Section IV presents the results of our 
calculations, together with a discussion of their appli¬ 
cability to the experimentally realizable compounds. A 
number of technical details are explained in appendices. 
Appendix A presents more details of the lattice struc¬ 
tures, their symmetries, and the band Hamiltonians of 
interest. Appendix B presents some technical details of 
the derivation of the Raman response. Appendix C gives 
detailed arguments for the polarization dependence of the 
Raman signal described in section IIIB. 


II. THE MODEL 

The pure Kitaev model that we work with below is a 
good starting point for understanding response functions 
in the spin liquid phase. The advantages of this model are 
threefold. First, it is exactly solvable on the whole family 
of harmonic honeycomb lattices. Second, both gapped 
and gapless spin liquid ground states can be obtained 
within this model by tuning the anisotropy of the spin- 
spin interactions on different bonds. Finally, the frac¬ 
tionalization of spin excitations emerges naturally in the 
exact solution,®® allowing for a clear identification of the 
role that these fractionalized excitations play in experi¬ 
mental response functions. For example, we show that 
only one type of fractionalized quasiparticle is probed by 
the Raman operator, which greatly simplifies the inter¬ 
pretation of the dynamical response. 

The Kitaev model®® on a generic tri-coordinated lat¬ 
tices takes the form 

n=Y. ( 1 ) 

(b>“ 

Here a = x,y and 2: label the three types of bonds eme- 
nating from each vertex, and and J° are the as¬ 

sociated coupling constants. 

Because only one component of the spin interacts along 
each bond, there is one conserved quantity for every 
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FIG. 1. (Color online) The harmonic honeycomb {'H-n) lattices: H-O (left), H-l (center), H-oo (right). Here n counts the 
number of rows along the c-axis before the orientation of the honeycomb plane switches between the two non-parallel chains of x 
and y bonds. The H-O Hyperhoneycomb lattice switches chains at every c-bond. The H-l has one set of c-bonds making rungs 
(gray) of ladders before a bridge (c-bonds in black) to the opposite ladder. The H-oo never switches ladders. The ladders are 
labeled by red(orange) and blue(green) on a;(t/)-bonds. The primitive unit cell for n < oo contains 4n-|-4 sites. The honeycomb 
lattice has only two sites in the unit cell. Our choice of unit cell is illustrated by the spheres, whose color alternates yellow and 
white indicating the odd and even sublattices respectively. The excited bonds in the tt flux state are circled on the H-l lattice. 


plaquette.^® This conserved quantity is given by 

Wp = l[af\ ( 2 ) 

ieP 

which is the product of spin operators around a plaquette 
P, whose spin component a(i) is given by the label of the 
outgoing bond direction, as illustrated for the honeycomb 
lattice in Fig. 2. 

The full spectrum of the model Eq. (1) can be de¬ 
scribed by using Majorana fermions to represent the 
spins. Following the notation of the original work by 
Kitaev,®° we use 

cr“ = iCjf^ . (3) 

The Majorana fermions cj and 6“ satisfy the algebra 

c 2=52 = 1, {6“,6f} = {6“,c,} = 0 (4) 

so that = c and Majorana operators at different sites 
anticommute. In terms of Majorana fermions the Hamil¬ 
tonian becomes 

H = i'^ (5) 

<b>“ 


where we have defined the bond operators 

= ibfb‘^ , ( 6 ) 

which satisfy Since the bond opera¬ 
tor commutes with the Hamiltonian in any eigen¬ 

state it is a constant. Therefore, in the unconstrained 
Hilbert space, eigenstates of Eq. (5) are described by the 
value of on each edge, together with an eigenstate 

of the resulting hopping Hamiltonian for the {cj} Ma¬ 
jorana fermions. Since = 1 then = ±1 

and the possible hopping parameters on each edge are 
tij — ^iJ■ 

Because the Majorana representation is redundant, not 
all of the eigenstates described above are physically dis¬ 
tinct. Imposing 1 = Dj = Cjbjhjbj at each site makes the 
representation exact. In particular, though individual 
commute with H, they do not commute with the 
constraint. The only conserved quantities corresponding 
to the in the physical Hilbert space are the {FFp}, 

defined by Eq. (2), which in the Majorana representation 
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FIG. 2. (Color online) An elementary plaquette of the hon¬ 
eycomb lattice with the bonds labeled by the component of 
spins that interact along them. The spin component shown at 
each vertex is the one that enters the product Wp of Eq. (2), 
and corresponds to the color of the bond pointing out of the 
plaquette from that vertex. 


in Fig. 1), a natural set of unit vectors (for n < oo) is®® 

ai = (-1,-^2, 0) 

a 2 = (-1,72,0) 

_ J (—1,0, 3)-I-(0, 0, 6 ) X I if n is even, , . 
\ (0,0, 6 ) X ^ if n is odd. 

where we have set the bond length to 1. The unit cell 
therefore contains 4n -b 4 sites. 

The spin-orbit coupling distinguished directions 
(x,y,z in Fig. 1), which determine which spin com¬ 
ponents couple along which bonds, are z = b, y = 
(a-b c)/-\/2, and x = (— a-bc)/-\/2. The nearest-neighbor 
bond vectors (in the a, b,c basis) are 


have the form: 

Wp= ■)<.). (7) 

{ij)ep 

Thus the physical eigenstates (and the band energies 
for {cj} Majorana fermions) depend only on the phase 
accumulated by a fermion hopping around a given pla¬ 
quette. We define the flux $p through plaquette P by 
exp[i4>p] = Y[^ij)^ptij/\tij\ = i'^Wp, where n is the 
number of the bonds around the plaquette P. Since 
= ±1, then Wp = ±1 and the fluxes - which ef¬ 
fectively generate lattice-scale magnetic flux for the dis¬ 
persing fermions - can only take on two values: 0 or tt. 
From here on we will therefore describe excitations of the 
Majorana operators in terms of their effect on these 
Z 2 fluxes, and refer to the dispersing Cj excitations as 
Majorana spinons. 


A. The Harmonic Honeycomb lattices 

Motivated by recent synthesis of the a-, fi- and 7 - 
Li 2 lr 03 compounds, here we study the Kitaev model on 
the "H-oo, 'H-O, and T-L-l lattices illustrated in Fig.l. The 
discussion here generalizes straightforwardly to the other 
lattices in the harmonic honeycomb series, j-L-n. 

The harmonic honeycomb series consists of bipartite 
orthorhombic tri-coordinated 3D lattices. The T-L-oo lat¬ 
tice is the exception, splitting into uncoupled 2D honey¬ 
comb lattice planes, shown on the right of Fig. 1. Be¬ 
cause the honeycomb lattice is well-known we leave the 
details of the in-plane coordinates and the unit cell that 
we use to Appendix A. The T-L-n lattice consists of n rows 
of co-planar hexagonal plaquettes, followed by a “bridge” 
layer of c axis bonds, and then another n rows in a new 
plane. In terms of the orthorhombic unit vectors (a, b, c 


d^ = ( 0 , 0 , 1 ) 

(black and gray) 

d^ = ^(l,72,-l) 

(red) 

d^ = ^(l,- 72 ,- 1 ) 

(blue) 

d^ = ^(-l,- 72 ,- 1 ) 

(orange) 

d^ = ^(-l, 72 ,- 1 ) 

(green), 


where the bonds are oriented from the odd (yellow) to 
the even (white) sublattice, and the colors refer to Fig. 1. 
Note that the x- and y-type couplings occur along bonds 
of two different orientations, labeled A (red x- and orange 
y-bonds) and B (blue x- and green y-bonds). All z bonds 
are parallel to the c axis. 

The spin-exchange Hamiltonian for each lattice is ex¬ 
pected to preserve the lattice symmetry, which, due to 
spin-orbit coupling, simultaneously acts on the spin and 
lattice basis vectors. On the 2D honeycomb lattice, Cg 
rotation symmetry^® requires that the exchange coupling 
be the same on each bond: . On the "H-n, 

n < 00 lattices the C 2 rotation symmetry interchanges J® 
and jy bonds, and thus . However, in general, 

jz ^ jx _ jy ^ since these are not related by symmetry. 
For comparison with the 3D cases, we will also consider 
jz ^ jx _ jy 2D lattice, although this explicitly 

breaks the Cg rotation symmetry. Moreover, the H-oo 
lattice does not have the full D^d point group symmetry 
of the 2D honeycomb lattice, because the vector relat¬ 
ing adjacent layers, ai, is not normal to the honeycomb 
plane (02 x c or ni x 712 (see Fig. 1)). This provides ad¬ 
ditional motivation for studying anisotropic couplings in 
these systems, even though in the limit that there is no 
interaction between the layers we expect that the full D^d 
symmetry will be restored.®^ A more detailed discussion 
of the lattice symmetries is given in Appendix Al. 

B. The spin liquid ground state 

For a given lattice, the ground state of the Hamiltonian 
(5) has the flux pattern that minimizes the energy of the 
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dispersing Majorana spinons. As finding the configura¬ 
tion of fluxes is lattice specific, here we will separately 
discuss the solution for each case. 

In the 2D honeycomb lattice {'H-oo) the ground state 
has zero flux.®° This can be proven exactly using a the¬ 
orem based on reflection positivity generally referred to 
as Lieb’s theorem. The theorem states that at least one 
ground state has flux <i>p = (n -I- 2)7r/4 through every 
n-sided plaquette that is symmetrically cut by a mirror 
plane that does not intersect any lattice points. It was 
originally proven for complex fermions in Ref. 71, see 
also Ref. 72. The extension to Majorana fermions was 
written down recently in Refs. 73 and 74. (In contrast 
to Ref. 73 we refer to the flux as the phase accumu¬ 
lated when hopping around a loop/plaquette following 
Ref. 50.) 

As explained in Appendix A, the 3D harmonic hon¬ 
eycomb lattices’ lack of mirror planes makes it difficult 
to uniquely fix the ground state flux configuration. Nu¬ 
merical searches®®’"^^ diagonalizing the Majorana spinon 
Hamiltonian in the thermodynamic limit for flux arrange¬ 
ments consistent with eight-fold enlarged unit cells®® or 
with 216 unit cells and periodic boundary conditions^® 
suggest that the ground state on the H-O lattice lies in 
zero-flux sector. However, on the 77-1 lattice a state with 
a simple, but non-zero flux pattern was found to have the 
lowest energy when = J^.®® We follow Ref. 69 

and call this the tt-Hux state. We have confirmed for 
the infinite lattice that this state is lower in energy near 
the isotropic point and that the zero-flux 

state becomes the ground state for small enough 
(we keep -P = P throughout) but that the energy differ¬ 
ence is very small (see Appendix A for the ground state 
energies). 

Once the flux pattern is fixed, the ground state is de¬ 
scribed by the band structure of the dispersing Majo¬ 
rana spinons in this flux background. Certain details of 
the band structure (though not the ground-state energy 
or other measurable quantities) are gauge-dependent, 
so to fix the band structure we work in a particular 
gauge. For the zero flux sector we choose the gauge 
“(d>“ “ ibfb'j = 1 for i on the odd sublattice and j 
even. If one represents the direction of positive by 

an arrow, then the arrows point from the odd-sites to 
even-sites.^® The other flux sectors are related by chang¬ 
ing the sign of the Majorana spinon hopping parameter 
on particular links. The tt flux state on the 77-1 
lattice can be obtained from this choice by changing the 
sign of W(y)“ on the bonds circled in Fig. 1. 

After finding the minimum energy flux sector and fix¬ 
ing an appropriate gauge, the Hamitonian (4) is reduced 
to a band Hamiltonian for the dispersing cj Majorana 
spinons. In the next subsection we discuss the important 
features of the resulting band structures, which can be 
probed directly with Raman scattering. 


C. The fixed-flux Hamiltonian 


We will now detail the fixed-fluxed eigenstates in 
terms of the band structures of the Majorana spinons. 
The Fourier basis is defined by Ck = 
where rg defines the position of the unit cell s, Cg = 
(cis, C3g,...; C2g, C4g...) defines the vector of distinct Ma¬ 
jorana spinons in this unit cell. Note that cj^ = and 
k belongs to the first Brillouin zone (BZ). 

In the presence of time-reversal symmetry, the Majo¬ 
rana spinon Hamiltonian on a bipartite lattice can be 
written in block off-diagonal form.®®d® We number the 
sites by their position in the positive c— direction with 
the odd-numbered sites (yellow) in the first block and 
the even numbered sites (white) in the second. Then 


H" = ^ i7k 


all k 


0 

-iK)' 


0 


= (ci.k, C3^k, C 2 ,k, C4,k, ••■)) (10) 


where 77 is a composite index 77 = n($) for the nP- 
harmonic honeycomb in the flux state 4) = {$p}p. If 
there are 2M sites in the unit cell, are M x M ma¬ 
trices. Their precise form for the cases we consider are 
given in Appendix A 4. 

To diagonalize the resulting band Hamiltonian, we ex¬ 
press the Majorana spinons in terms of complex fermion 
operators as 


Ck — V^Qk^^k 

Q^k “ (Ol.ki fl2.k) j ^i.-kJ ®2,-k> •■•) “ (^ki^-k)’ (H) 

where the o„.k satisfy the usual complex fermion algebra 
k’ ®m,q} = bnmSkq and Qk is a unitary matrix that 
diagonalizes the matrix H^: 

H" = ^ a^k^k«k = 2QIhIQ^. (12) 

all k 

The diagonalized Hamiltonian reads as 

= 2 'y ( ~ Qm,koJ„_k)’ (1^) 


where m labels the bands. 

The fact that j, is related to a„,k by hermitian con¬ 
jugation leads to some redundancy. This, along with 
the condition cj^ = c_k, requires that = yQlk’ 

where 7 = f ? n V In practice we impose this con¬ 


straint by computing Qk for k in half of the Brillouin 
zone and defining (5_k = iQ'L- The ground state is the 
one with no quasi-particles, am,k |0) = 0, and has energy 

For every two Majorana spinons there is one spin¬ 
less fermion, whose corresponding excited state |0) 
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has energy £m,k- We refer to the energies of the differ¬ 
ent quasiparticles as different Majorana spinon bands, of 
which there are one for every two sites in the unit cell. 
The physical excitations above the ground state corre¬ 
spond to pairs of the original Majoranas: flux-|-Majorana 
spinon or two Majorana spinons. Operators that cre¬ 
ate single-Majorana fermion excitations do not commute 
with the constraint, and therefore are not physical. 

The band structures that result from this picture have 
some striking similarities for all three of the lattices 
considered here. First, there are two distinct spin liquid 
phases: a gapless phase centered at the isotropic point 
jx — jy — jz^ gapped phase when > (J^^-l-J^). 

Second, the codimension of the Fermi surface in the 
gapless phase is always 2:^®d8.79 jg ^ Dirac line in 3D, 
and a Dirac point on the 2D honeycomb lattice. As a 
result, in the gapless phase, the low-energy density of 
states asymptotically obeys the power law p{ui) ~ w. 

III. RAMAN RESPONSE 

In this section we present the derivation of the Raman 
response and the polarization dependence of the Raman 
spectrum for the Kitaev spin liquid state on three har¬ 
monic honeycomb lattices: the "H-oo, "H-O, and %-! cor¬ 
responding to a-, /3-, and 7-Li2lr03, respectively. While 
our study does not represent a complete microscopic the¬ 
ory of Raman spectra (including all possible contribu¬ 
tions), it should provide experimentalists with the main 
signatures of magnetic quantum number fractionalization 
that might be visible in these compounds. 

Raman scattering is the inelastic scattering of light by 
energies in the meV range. It measures correlations be¬ 
tween two-photon events (one ingoing and one outgoing). 
By Fermi’s golden rule the Raman intensity can be writ¬ 
ten as the Fourier transform of a correlation function 

I{w) = J (i?(f)R(0)), (14) 

where R{t) = is the Raman operator in the 

Heisenberg representation. We derive the Raman opera¬ 
tor using the Loudon and Fleury (LF) approach for Mott 
insulators,which describes the effective interaction of 
light with spin degrees of freedom. The dominant cou¬ 
pling was shown to be the electric one.®° In these sys¬ 
tems, the effective interaction can be obtained by per¬ 
forming a hopping expansion; the leading term is the LF 
Hamiltonian.^® 

If the incoming photon frequency is smaller than the 
appropriate Mott gap,®^ the resulting LF Hamiltonian 
turns out to be precisely the exchange Hamiltonian aug¬ 
mented with polarization-dependent terms correspond¬ 
ing to the component of each photon’s electric field along 
the bond that the electron virtually hops across. For a 


spin-exchange Hamiltonian H this leads to®®’®^ 

(15) 

R= ^ (ei„.d,, )(£;„. d,,)F“^u“af (16) 

where ein and Cout are the polarization vectors of the 
incoming and outgoing light, is the vector from site i 
to site j, af is the component of the spin on site I, 
and F is a generalized exchange constant. 

A. Quasiparticle picture 

As emphasized in the introduction, one striking feature 
of Raman spectroscopy in a Kitaev spin liquid is that the 
Raman operator couples only to the dispersing Majorana 
spinons. For the Kitaev Hamiltoian (1), the LF operator 
Eq. (16) takes the form®^ 

R=Y. (^in • d“)(eout • d“)J“a“a“ 

(b>“ 

= i ^ (cin • d“)(eout • d“) = (17) 

(b>“ 

which is a simple quadratic operator in terms of the Ma¬ 
jorana spinons Ck. Remarkably, due to its similarity to 
the original Kitaev spin Hamiltonian, the Raman opera¬ 
tor (17) does not excite the gapped flux excitations. This 
is one distinct advantage of using Raman response in 
these systems, as it probes only one of the fractionalized 
sectors. (In contrast, neutron scattering always excites 
both fluxes and Majorana spinons.) In addition, the fact 
that the Raman operator conserves the flux in each pla- 
quette greatly simplifies the calculations of the Raman 
intensity, since we can use the fixed-flux Hamiltonian 
(10). Consequently, we can write the Raman operator 
in terms of the Bogoliubov-deGennes fermions that diag¬ 
onalize the Hamiltonian (see Eq. (11)). In this basis, the 
Raman operator takes the form R = i?k with 

R\^ — 4“ 2 —k 4“ h.C.^ , 

(18) 

where indices m,n = 1,...,M refer to bands and Amn,k 
and Bmn,k are bilinear functions of the in and out po¬ 
larizations. One can reduce the Raman intensity to a 
weighted two-Majorana spinon DOS (2-DOS) given by 

P2{0J) = ^ d(w - £m,k - En.k)- (19) 
Specifically, we obtain 

/(uj) = TT ^ ^ £'m,k ^n,k) |-^mri.,k| 

m,n;k 



7 


Note that is the antisymmetric matrix correspond¬ 

ing to the creation of two excitations in Eq. (18). The 
^mn.k part annihilates the ground state and thus does 
not contribute to the Raman response. 

Eq. (20) describes the response at zero temperature. 
At finite temperature the Amm',k can contribute to the 
correlation function by relaxing thermally excited quasi¬ 
particles. In general the finite temperature dependence is 
complicated by the presence of thermally excited fluxes, 
which are strongly coupled to the Majorana spinons. 
However, for temperatures well below the flux gap one 
can restrict the calculation to the zero flux sector. In 
this case Eq. (20) generalizes to 

/(cu) = 27r ^ ( s kTlnmjk*^ (^ “b [^m,k ^n,k]) 

m,n;k ^ 

X nF(£m,k)n_F(-e„,k) 

”^2 |-^mn,k| ^ (^ [£m,k £n,k]) ^f( £m,k)^F( £n,k)^ 5 

( 21 ) 

where npis) = is the Fermi-Dirac distribution. 

B. Polarization dependence 

In the following we will identify the distinctive features 
of the polarization dependence of the magnetic excita¬ 
tions in the Kitaev model and their relation to the lat¬ 
tice structure and the form of the electronic Hamiltonian. 
But first we review the arguments for the polarization 
dependence that follow directly from lattice symmetry 
and the form of the LF operator. Due to the general¬ 
ity of these assumptions, the relationships obtained are 
expected to hold for any nearest-neighbor spin-exchange 
Hamiltonian on these lattices, such as the Heisenberg- 
Kitaev model. Ways to break these constraints on polar¬ 
ization dependence are discussed in the next section. 

The form of the Raman operator R depends on the 
polarization of incoming and outgoing light. For the LF 
operator [Eq. (16)] this is because driving an exchange 
process along a certain bond requires the photon’s polar¬ 
ization to have a component along that bond vector. We 
will use the short hand notation /if to refer to a scat¬ 
tering geometry with the polarization of incoming light 
along /t and outgoing along v. It is convenient to use the 
cubic coordinates /i, f = a, &, c (see Fig. 1) due to their 
relation to the symmetries of the lattice (see Appendix 

Al). 

Since the Raman operator is linear in polarizations it 
can always be written as a tensor dotted with the polar¬ 
ization vectors. 


are 

(23) 

a=bonds 

where a refers to a particular bond direction, (i“ is the 
component of the bond vector (i“, and H°‘ is the sum 
of the spin-exchange terms on bonds directed along a. 

By inserting Eq. (22) into Eq. (14) the Raman intensity 
can be decomposed into a linear combination of intensi¬ 
ties obtained for pairs of these Raman operators. 

J(w) = J {R{t)R{0)) 

= ^ ^ (^in)/j-(^out)i/(^in)/i'(^out) (^) 7 

,iy' 

(24) 

where 

V,^,,,(w) = J ((R^.(t)RF'-'(O)) 

+ {R^>,'{t)R^,{{)))) (25) 

The terms with fJ,,v d' i^' only appear when either the 
incoming or outgoing polarization is not along a cubic 
unit vector a, &, or c so that the Raman operator is the 
sum of multiple terms. For instance, a polarization with 
Cjn = a and Cout = ^ would give 

d- — 2 {daa^aa “b q^) \/^Iaa,ab: (26) 

where we have left the to dependence implicit, to sim¬ 
plify the expression. One could measure the spectrum 
Iaa,ab by Hrst isolating the other terms Iaa,aa and Iab,ab 
and then subtracting off the components of these other 
two from the spectrum measured with the polarization 
in Eq. (26). 

Note that the intensity “cross terms” Iab,cd, ab ^ cd are 
sums and differences of measurable quantities and may 
therefore be positive or negative. Otherwise, we refer 
to the intensity for a simple polarization configuration 
along the cubic directions fj, and v by which is given 
by above. 

The space group symmetries reduce the number of 
independent spectra, which can be classified by the ir¬ 
reducible representation (irrep) under which the point 
group acts on the Raman tensor The cross terms 

between different irreps are zero since symmetry must act 
trivially on the spectrum I. As discussed in more detail 
in Appendix C, the point group D 2 h reduces the number 
of non-zero spectra of the 3D lattices to nine (because 12 
cross terms vanish): 


R — 'Y, i^in)yiRyLiyi^ont)v ■ ( 22 ) 

For a spin-exchange Hamiltonian the Raman operators 


^bb-i ^cct ^aa^bbi ^aa,cci ^bb^co ^abi ^bc 7 ^ 0 - 


(27) 


The first three correspond to different representations of 
the Aig irrep, the next three are their cross terms, and 



the remaining three correspond to the Big, i? 2 g and B^g 
irreps respectively. 

In addition to the point group symmetries, the 3D har¬ 
monic honeycomb lattices have an additional constraint 
that can be understood in terms of a non-symmorphic 
screw symmetry of a related lattice of the same connec¬ 
tivity. This reduces the number of independent spectra 
to 6 by the relations: 

~ ^bb-} ‘^^aa,cc ~ ^bb,cci ‘^^ac ~ ^bc- 

We discuss this in detail in the results section below. 

Quite separately, the close resemblance of the LF oper¬ 
ator Eq. (16) to the spin-exchange Hamiltonian can give 
additional linear relations between the spectra. This is 
due to the fact that the Hamiltonian acts trivially on 
eigenstates. Therefore if two Raman operators sum to 
the Hamiltonian, their spectra must be related. We refer 
to this as the Loudon-Fleury (LF) relationship. For the 
remaining six spectra, this type of argument leads to 

Icc = QIaa = ~3/oo,cc (29) 

for a nearest-neighbor spin-exchange Hamiltonian on the 
3D lattices. That leaves only three independent spectra 
for the model considered here. We choose to plot laa, lab, 
and lac, which represent the Aig, Big, and B 2 g channels. 
Then the other channels are either zero (see Eq. 27), or 
are linearly related to these ones by Eqs. (28-29). In 
particular, B^g is degenerate with i? 2 g- 

Eor the 2D lattice with full lattice symmetry there 
are two independent spectra by point group symmetry 
and one LF-relation {RAig = H). This reduces it to 
one independent spectrum - the Eg channel, which we 
represent with Ixx-^^ In Appendix C we review this ar¬ 
gument as well as treating the case of bond anisotropy 
jz ^ jx _ jy breaks the Cq rotation symmetry. 

We find that there are then four independent spectra by 
symmetry and the same LF-relation reduces this to two. 
In that case we plot Ixx and Ixy which represent the Aig 
and the Big channels respectively, of D 2 h- Inclusion of 
next-nearest-neighbor spin-exchange couplings does not 
change the vanishing of the Aig channel due to the LF- 
relationship, although it would break the two equalities 
in Eq. (29) for the 3D lattice. See Appendix C for more 
details. 


IV. RESULTS 

In this section, we present our results and discuss the 
key features of the Raman response of the Kitaev model 
on the H-O, H-l, and "H-oo lattices. For all three lat¬ 
tices the Raman response differs qualitatively from the 
structure factor. This follows from the fact that Raman 
directly probes the fermion DOS, without coupling to 
fluxes. The difference is especially pronounced in the 
gapless phase, where the Raman spectrum is gapless but 


the structure factor is gapped. In addition, we will high¬ 
light certain distinctive features of the polarization de¬ 
pendence in the Raman spectra for a Kitaev spin liquid 
phase on these lattices. 


A. 2D Honeycomb spectra 

The simplest case is the 2D honeycomb lattice. It has 
one Majorana spinon band and at most four independent 
polarization combinations by symmetry. At the isotropic 
point = jy = the Cq rotational symmetry allows 
only two independent Raman spectra corresponding to 
the channels Aig and Eg, which are probed by the com¬ 
binations 

^Aig = 2 i.^xx ~ Ixy) (30) 

lEg = Ixy (31) 

As discussed in the previous section, there is no response 
in the Aig symmetry channel (it does not couple to the 
spins) due to the LF relationship. For equal couplings, 
the rotation symmetry gives a polarization independent 
spectrum corresponding to the Eg channel. This is a pe¬ 
culiar prediction that holds for any symmetry-preserving 
spin-exchange Hamiltonian with a ground state that also 
preserves the lattice symmetry. In particular it expected 
to hold in the Kitaev spin-liquid phase in the presence of 
a Heisenberg exchange perturbation. The resulting spec¬ 
trum for the pure Kitaev model is plotted in Fig. 3(a). 

When we take = J^, the rotation symmetry is 

broken and there are four independent non-zero spectra. 
Nevertheless, they are reduced to two because the LF 
relationship forces the three terms in the Ag-channel to 
be degenerate {Ixx = lyy = —Ixx,yy)- The two remaining 
spectra are represented by Ixx and Ixy of the Aig and Big 
irreps. These are plotted for a gapless point {J^ = = 

1.43 J^) and a gapped point (J“ = = 0.3J^) in Figures 

3(b) and 3(c). The breaking of the LF relationship is 
measured by a nonzero sum Ixx + Ixx,yy in all cases. 

Along with the excitation gap seen in the gapped 
phase, the computed intensities show broad humps that 
are qualitatively similar to the isotropic case. The rela¬ 
tive total spectral weight between the two active chan¬ 
nels is plotted in Fig. 4 where we find that the rel¬ 
ative weight is determined by the coupling anisotropy 
. One can understand this dependence by consider¬ 
ing the following limits: (1) = 0, where the states 

are zero-dimensional, with no Raman response; and (2) 
-A oo, where the states are one-dimensional, with 
Raman intensity only in the A 2 g channel. 

Apart from the broad hump there are also fine fea¬ 
tures related to the Majorana spinon density of states 
(shown in the left insets). In particular, the van Hove 
singularities of the saddle points in the dispersion lead 
to singularities in the derivative of the Raman response 
(shown in the right insets of Fig. 3). 
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FIG. 3. (Color online) Raman intensities for the pure Kitaev model on the 2D honeycomb lattice computed for: (a) the 
symmetric point = J^, (b) another gapless point = 1.43J'', (c) a gapped point = 0.3J^. For each 

figure, the inset on the left is a plot of the one-particle DOS; the inset on the right is a plot of the derivatives of the Raman 
intensity I' {oj). 



FIG. 4. The relative total spectral weight SW = J I{Lo)du! 
of the two representative Raman intensities as a function of 
the Kitaev exchange-coupling anisotropy for the 2D honey¬ 
comb lattice. The transition to the gapless phase occurs at 
r/r = 0.5. 


B. 3D Harmonic Honeycomb spectra 

Many of the characteristic features of the 2D Raman 
spectra carry over to the 3D lattices. Notably, because 
Raman scattering couples only to the Majorana spinons, 
the spectra are broad, reflecting the total bandwidth®^ 
AzJ°‘^ = 4(J® -\- jy -\- J^) accessible to two-particle Ma¬ 
jorana spinon excitations above the ground state, where 
z is the number of neighboring sites. We note that this 
is also the energy range in which we would find magnons 
(spin-waves) if the state were ordered. Moreover, in the 
gapless spin liquid phases we always find a linear Raman 
spectrum at low energies, as can be anticipated from the 
ID nature of the fermi-surface (a line) and the linear dis¬ 
persion away from that line.®®’^® This feature, together 
with a gapped dynamical spin structure factor, are com¬ 
mon to all of the harmonic honeycomb lattices. 

However, there are also a number of notable differ¬ 
ences. The most striking one is the greater number of 
independent, non-vanishing spectra as a function of po¬ 


larization, due to the lower symmetry and larger number 
of ways for light to couple to the lattice. This can be 
clearly seen in the Raman spectrum of the "H-O lattice, 
shown in Fig. 5. In particular, we find that 2Iac = Ibc so 
that the i? 2 g and B^g channels have identical response. 
Although no space group symmetries relate the a and b 
directions, this result can be understood in terms of the 
symmetry of an equivalent lattice model. The connec¬ 
tivity is unchanged by scaling a —>■ v^a, which leads to 
a lattice with screw symmetry composed of C 4 rotation 
about an inversion center, along with 03/2 translation. 
Equivalently, the original lattices are symmetric under 
the scaling a V^a,b —>■ b/V^ followed by the same 
screw rotation. This symmetry exchanges a and b as 
desired, without affecting the physics of the spin model. 
(We use this transformation only as a theoretical tool: in 
practice, macroscopically compressing the Li 2 lr 03 mate¬ 
rials would break the symmetry between the spin-orbit 
distinguished directions - see Refs. 57 and 68 for material 
details). The factor of 2 between the intensities lac and 
Ibc comes from the different projections along the actual 
a and b directions. In addition, the equalities Ilaa = hb 
and 2 Iaa,cc = Ibb,cc mentioned in the previous section are 
also guaranteed by this symmetry, leaving only two rela¬ 
tionships that depend on the Loudon-Fleury form of the 
operator: I^c = 9Iaa = -3/aa.cc- 

Unlike in the 2D case, the low-frequency response for 
the H-O lattice shows a strong polarization dependence, 
with only one of the three polarization combinations 
showing a significant low-frequency response in the gap¬ 
less phase. In fact, a strong polarization dependence in 
the low-energy response is common to all of the 3D lat¬ 
tices H-n, n < 00 : the Big channel {lab in the Figure) 
is always inactive at low energy. This can be understood 
by considering the symmetry of the Raman operator in 
this channel. There is a Z 2 glide plane symmetry (see 
Appendix A I) common to all finite n lattices. Therefore 
the bands can be labeled even or odd under this sym¬ 
metry. For momentum-dependent symmetries such as a 
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FIG. 5. (Color online) Raman intensities for the H-O lattice computed for: (a) (gapless), (b) — 1.43J^ 

(gapless), (c) — 0.3J^ (gapped). The spectrum looks qualitatively the same throughout the gapless phase, in contrast 

with the 2D case in which the symmetric coupling point has extra symmetry. The maxima and minima in I{lu) originate 
from van Hove singularities and band edges, respectively, for each band (see Fig. 6). The only two particle combination that 
couples to the Big channel is one from the lower band and one from the upper band causing the spectrum lab to vanish unless 
4J" < a; < 4^5 JC®'* 



FIG. 6. (Color online) The Majorana spinon spectrum on the B-O lattice plotted along high symmetry lines. The inset on the 
left shows the 2-DOS defined in Eq. (19). For details of the BZ see Appendix A3. The two complex fermions per unit cell make 
two bands, e±,k. This gives three distinct configurations for the two-particle states of interest in Raman scattering: two with 
both excitations in the same band, and one with one excitation per band. The 2-DOS is split into the contributions from each 
of the configurations by pm.n{^) = “ £m,k — £n,k), where indices m,n = ±. The maxima and minima in /(w) appear 

as extrema in the spectrum of the appropriate band at some symmetry-distinguished point in the BZ. 


glide plane, the parity of the Raman operators can be 
carried by momentum-dependent coefficients rather than 
the excitations. However, the Raman operator in the Big 
channel {Rab) is odd under the glide plane transforma¬ 
tion, and, since this operator does not vanish at k = 0, 
this parity must be carried by the Majorana spinon exci¬ 
tations that it creates. Therefore the pair of excitations 
must change the glide plane symmetry quantum number 
of the state it acts on. Because glide plane symmetry 
is a Z 2 symmetry, creating two Majorana spinon excita¬ 
tions in the same band is always a glide-even operation. 
Therefore the Raman operator in the Big channel can 
only couple to two-particle states involving excitations 
from different bands. In particular it cannot excite a pair 
of excitations both in the lowest band. (Nor can it ex¬ 
cite a pair in the highest band). This is apparent in Fig. 
5(a) where Rb (the Raman intensity in the Big channel) 
is non-vanishing only in the region 4 < to < - 

i.e. only at energies accessible by exciting one Majorana 
spinon from each band.®^ For identification of the contri¬ 
butions of each band see Fig. 6 below. Since for general 
n the Majorana spinon ground state contains only one 
gapless band, the spectrum lab is gapped for all of the 
3D lattices. See Figures 7 and 8 for similar vanishing 


of lab on the H-l lattice. This is a particular result of 
the interplay between the fermionic fractionalization in 
the Kitaev spin liquid and the symmetries of the 3D har¬ 
monic honeycomb lattices. 

Another notable difference between the Raman spec¬ 
tra of the 77-0 and R-oo lattices is the number of maxima 
and minima. The 2D lattice shows only a broad contin¬ 
uum (at least in the pure Kitaev model®''^), with only one 
broad peak. In contrast, the 3D lattices show a number 
of rather sharp peaks in addition to the overall broad 
response. This is even more noticeable in the spectrum 
of the 77-1 lattice, shown for the zero flux ground state 
in Fig. 7, and for the tt flux ground state, which has an 
enlarged unit cell, in Fig. 8, for = jy = . 

These new peaks can be understood by considering 
the momentum-locked two-particle density of states (2- 
DOS) defined in Eq. (19), which the Raman spectrum 
closely emulates. For the 77-0 lattice, this is plotted in 
Fig. 6, along with the band structure along high symme¬ 
try lines. The 4-site Majorana spinon unit cell implies 
that there are two Majorana spinon bands, and there¬ 
fore three types of quasiparticle pairs above the ground 
state. As the figure makes clear, the sharp features cor¬ 
respond either to band edges, or to van Hove (VH) sin- 
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gularities that occur near the high-symmetry points in 
the Brillouin-zone. 



FIG. 7. (Color online) Raman intensities for the H-l lattice 
in the zero-flux state at the isotropic point {J^ = = J^). 



FIG. 8. (Color online) Raman intensities for the H-l lattice 
TT-flux state at the isotropic point (J^ = = J^).®® 

Fig. 9 compares the 2-DOS of the H-O, H-l, and H-Itt 
lattices. (Here H-ln refers to the tt flux ground state 
of the H-l lattice, while H-l denotes the energetically 
proximate 0-flux state.) Since the number of sites in the 
unit cell increases with n, so does the number of sharp 
features in the Raman spectrum (see Figs. 3a, 5, 7, and 
8 ). In the 2D lattice there is only one Majorana spinon 
band, and the Raman spectrum has only an upper edge, 
a lower edge and a VH singularity. For the H-n lat¬ 
tice there are n -|- 2 bands. Then each pair of Majorana 
spinons has a maximum and a minimum, as well as possi¬ 
bly one or more VH singularities, leading to a number of 
features growing approximately quadratically in the unit 
cell size. 

To summarize, the Raman spectra of the 3D harmonic 
honeycomb Kitaev models differ from those of the 2D 
model in both their polarization dependence and the 
number of sharp peaks. However, a number of key fea¬ 
tures are common to all of these spin liquid phases, 
including the bandwidth and the fact that symmetry 
and LF relations determine the polarization dependence. 
Moreover, in the 3D lattices the polarization dependence 
is again characteristic of the coupling anisotropy 


This is illustrated in Fig. 10, which shows the relative to¬ 
tal spectral weights (integrated over all frequency) of the 
representative polarization combinations for the H-0 lat¬ 
tice. We have not plotted the analogous weights for the 
H-l lattice because the results are nearly indistinguish¬ 
able. The limits of j = 0, oo can be understood 
the same way as for the 2 lattice. For the 3D lattices 
we do not have polarization independent response at the 
symmetric coupling point as there is no 

change in symmetry at that point. 


C. Discussion 

Our discussion so far has focused on the exactly solv¬ 
able Kitaev Hamiltonian at zero temperature. While this 
has the advantage of allowing us to evaluate the Raman 
spectrum exactly, in any real experiment additional spin- 
spin couplings would be present, which would have to be 
taken into account to predict the Raman response even 
within the spin liquid phase. In addition, it is useful 
to understand how finite temperature affects the spectra 
shown above. Here we will comment on how our results 
are expected to change at finite temperatures, and in 
the presence of perturbations away from the pure Kitaev 
Hamiltonian. Our discussion is qualitative rather than 
quantitative in nature, due to the technical challenges 
involved in repeating our calculation in these regimes. 

First, we consider how the presence of additional spin- 
spin couplings would affect our result. In general, such 
interactions will not commute with the conserved quan¬ 
tity Wp (Eq. (2)). As a result there will be some flux 
fluctuations in the ground state, which is no longer ex¬ 
actly described by a Majorana spinon band structure. 
However, the most noticeable change is likely the one due 
to the change in the LF operator, which will have some 
finite probability to excite fluxes. Ref. 51 showed, us¬ 
ing a perturbative treatment of a small nearest-neighbor 
Heisenberg interaction, that this leads to an additional 
peak in the Raman response at the energy to create four 
flux excitations. For the 3D lattices we expect the situ¬ 
ation to be qualitatively similar: adding a small Heisen¬ 
berg (or other symmetry-preserving) perturbation to the 
Kitaev Hamiltonian should lead to a new peak at the en¬ 
ergy scale of the gap for two flux-loop excitations. We 
note that quantitatively this flux gap for the 3D lattices 
is expected to be smaller with respect to the Kitaev ex¬ 
change coupling than in the 2D case.^® 

As discussed in Section HIB, the polarization depen¬ 
dence is derived from two considerations. Many of the 
relationships between polarizations are fixed by lattice 
symmetry alone. However, some of the relationships be¬ 
tween different polarization channels found in Appendix 
C come from the close relationship between the Hamilto¬ 
nian and the LF Raman operator, which we call the LF 
relationship. In doing perturbation theory, it is there¬ 
fore essential to treat the Hamiltonian and the Raman 
operator to the same order to preserve this relationship. 
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FIG. 9. 2-DOS for the Majorana spinons on the (a) T-i-Q, (b) H-l and (c) H-ln lattices. Up to polarization dependence these 
spectra determine the character of the Raman intensities. The number of maxima and minima grows linearly with the number 
of combinations of bands for two-particle states, which increases from left to right due to increasing unit cell size. 



FIG. 10. (Golor online) Relative spectral weights for Ra¬ 
man on the H-O lattice as a function of the bond-coupling 
anisotropy. 


We note that Ref. 51 treated corrections to the Raman 
operator but not to the Hamiltonian (and therefore the 
ground state) in the presence of a Heisenberg pertur¬ 
bation. They consequently found that this term breaks 
the polarization independence described above. How¬ 
ever, adding perturbative corrections to the same order in 
the ground state should restore this symmetry, and pro¬ 
duce the polarization-independent Raman response an¬ 
ticipated from the combination of symmetry and the LF 
relationship. 

The LF relationships reported here are expected to 
hold for all nearest-neighbor spin-exchange Hamiltoni¬ 
ans. Breaking an LF relationship results in more linearly- 
independent spectra. On the 2D lattice the same LF re¬ 
lations apply for any bond orientation. However, for the 
3D lattices the LF relationships can be broken by the in¬ 
troduction of further neighbor spin-exchange terms in the 
Hamiltonian. See Appendix C for more details. There 
are also ways to break these relations that apply to all the 
lattices. First, perturbations that do not come from an 
electron virtual hopping process - such as a time-reversal 
(TR) breaking magnetic field - would break this relation¬ 
ship. Also, tuning the incoming photon energy into the 
resonant regime (near the Mott gap) will allow multiple 
electron hoppings and therefore leads to the appearance 
of higher order spin terms in the Raman operator that 


are not present in the effective Hamiltonian.^®’®^ This is 
also expected to break the LF relationship. We also note 
that a small TR-breaking perturbation can also lead to 
a change in the power law of the low energy Majorana 
spinon DOS.^® We will discuss this in the context of Ra¬ 
man scattering in future work.®® 

Next, let us consider the effects of finite temperature. 
At temperatures much lower than the flux gap, flux ex¬ 
citations are rare, and the spectrum is expected to look 
very similar to the zero-temperature result derived here. 
The transition out of the spin-liquid phase occurs when 
flux loops proliferate, confining the dispersing Majorana 
spinons. The transition is expected to be different in 
the 2D and 3D cases, due to the different nature of the 
flux excitations. In 2D the flux defects are point-like, 
with no long-ranged flux-binding interactions. At any fi¬ 
nite temperature there will therefore be a finite density 
of unbound fluxes, which immediately confines the Ma¬ 
jorana spinons at long length scales due to the mutual 
semionic statistics between Majorana spinons and fluxes. 
In 3D, however, the flux defects are loops, with an energy 
cost proportional to the total loop length E oc aL, where 
length L is measured in plaquettes. Below Tc, which is 
set by a,^® the flux loops stay small, and since in 3D the 
mutual statistics is only felt by Majorana spinons when 
they go through a loop, a finite density of small loops 
does not lead to confinement. At temperatures above Tc 
flux loops may still be rare but the entropy gain in hav¬ 
ing larger loops wins over the energy cost and there is a 
small density of loops of arbitrary length.^®’®® 

Above the critical temperature there is no spin liq¬ 
uid phase because the Majorana spinons are confined de¬ 
grees of freedom and therefore cannot form bands. How¬ 
ever, there is a confinement length scale (set by the den¬ 
sity of loops) below which the Majorana spinons can still 
propagate. This sets an energy scale above which the 
Majorana spinons may have dynamics that reflect the 
deconfined (spin liquid) state. The finite temperature 
transition has been considered in quantum Monte-Carlo 
for the H-O lattice.®®“®® The transition temperature turns 
out to be about Tc = O.OIITr- for the H-O lattice with 
jx _ jy _ jz ^ where Tk is the temperature correspond¬ 
ing to the Kitaev exchange coupling J“. In fact, Ref. 
88 claims that the dominant changes appear only in the 












13 


low-energy Majorana spinon DOS as the temperature in¬ 
creases in the range Tc < T < Tk- Therefore, the high 
energy features of our Raman response could be observ¬ 
able up to < Tk, which is also consistent with recent 
experiments.®^ 

Finally, let us speculate on the possibility of detect¬ 
ing spin-liquid signatures in the Raman spectra of the 
A 2 lr 03 harmonic honeycomb irridates, and a-RuCls. 
Experimentally, these materials have been found to order 
at low temperatures.®® This is not unexpected, since the 
appropriate spin Hamiltonian contains both Kitaev inter¬ 
actions and other spin-exchange terms.However, 
one might hope that signatures of a proximate spin-liquid 
phase are observable at temperatures above the ordering 
temperature Tk- Specifically, if the quantum transition 
out of the spin liquid phase due to perturbations is first 
order, finite temperatures fluctuations for T > Tk will 
carry signatures of the spin liquid state. This case has 
already been made in Ref. 92 - see also Ref. 3. 

It is therefore not unreasonable to expect that the Ra¬ 
man spectra of the currently available compounds in the 
temperature range Tk < T < Tk could display high- 
energy features of the Kitaev spin-liquid state. For in¬ 
stance, after isolating the magnetic contributions from 
the ones due to phonons, one could look for the band- 
resolved polarization dependence found for the 3D lat¬ 
tices above. As already mentioned, much of the polariza¬ 
tion dependence is dictated by symmetry for any spin- 
exchange Hamiltonian whose eigenstates respect the lat¬ 
tice symmetry. Overall, we believe that any quantitative 
agreement with the high-energy Raman spectra reported 
here as well as with the distinct polarization dependence 
indicative of QSL phases on the harmonic honeycomb 
lattices would make the case for these long sought-after 
states of matter much stronger. 

In summary, we have studied Raman scattering in the 
Kitaev spin model on the harmonic honeycomb lattices. 
We hnd that the Raman response of the spin liquid has 
the following distinctive features: 

1. A broad continuum response determined by the 
Majorana spinon 2-DOS and specific features char¬ 
acteristic of the unit cell size. 

2. A gapless phase with a low energy asymptotic 
power law reflecting the exotic Fermi surfaces. 

3. Linear relations between distinct polarization chan¬ 
nels dictated by the LF relations. 

4. Band-resolved polarization dependence in certain 
channels for the 3D lattices that reflects the partic¬ 
ular fractionalization of the Kitaev spin liquid. 

5. Polarization dependence of the spectral weights 
that is related to the exchange-coupling anisotropy. 
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Appendix A: Harmonic Honeycomb lattices 

Here we provide additional details on the structures, 
symmetries, and Brillouin zones of the harmonic honey¬ 
comb lattices. 

We Hrst introduce some vocabulary. On each lattice 
the X and y bonds make chains along two distinct direc¬ 
tions labeled A and B in Eq. (9). We call the bonds that 
connect parallel chains rungs and the ones connecting op¬ 
posite chains bridges - see Fig. I. The 2D honeycomb 
lattice has only rungs making up a single ladder. In the 
c-direction, the B-n ladders have n -I- I chains connected 
by n sets of rungs before a bridge to the opposite ladder. 

For the 3D lattices, the coordinates that we use are 
given in the main text. (See Eqs. (8) and (9)). For the 
B-oo lattice (2D honeycomb) it is more convenient to use 
coordinates in the lattice plane as in Ref. 51. For this 
case we rename c as y and define x = B oc (a, 6,0). We 
reserve square brackets for vectors in these coordinates. 
The unit vectors are ni /2 = [3,±-\/3]/2 and the bond 
vectors (from odd sublattice to even) are given by 


= [0,1] 

(blue) 


rf" = ^[-V3,-l] 

(red) 



(green). 

(Al) 


1. Lattice symmetries 

Here we review the symmetries of the three lattices, 
which are used in Appendix C to constrain the polariza¬ 
tion dependence, and also below to constrain the ground 
state flux conhguration. 

The n < oo B-n lattices all have G 2 symmetries in 
the orthogonal a, b, and c directions about the center- 
points of bridges. Note that the a and c ones swap x 
and y bonds, so that if ^ the point group would 
be broken down to C 2 h, although we do not consider 
that case here. They also have inversion centers at the 
midpoint between the center-points of two near bridges. 
All of them have glide planes (reflection x translation by 
half a unit vector) in all three of the directions a, b, c. For 
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the odd-n lattices, the glide planes are a, 6 , or c-reflection 
about a bridge-center with 03/2 translation. The odd-n 
lattices also have mirror symmetries in the c-direction 
about the inversion centers. 

For the even-n lattices the glide planes’ reflections pass 
through the inversion centers, which sit on the center of 
an a; or y bond at the midpoint between two bridge cen¬ 
ters. The glide plane c —)■ —c must be composed with 
a 2 -translation or oi-translation for centers of A- or B- 
type bonds respectively. The glide plane about a on 
and bonds is composed with (03 -I- ai)/2 translation 
while on J\ and bonds it includes (03 -I- a 2)/2 trans¬ 
lation. Finally, the glide plane taking b —>■ —b includes 
03 /2-translation or (oi -I -02 + 03)72 translation when the 
reflection plane passes through x- or j/-type bonds respec¬ 
tively. 

The 2D honeycomb lattice is more constrained by sym¬ 
metries, including Cq rotation and mirror planes.™ The 
3D point group of the 2D lattice is Dsd- The projection 
of the point group onto the 2D plane is generated by 
{Cg, cr} where 




1 

V3 1 ) 



(A2) 


In coordinates {ka,kb,kc), the high-symmetry points 
used in Fig. 6 are 


F= (0,0,0) 

/ IItt tt 

— 1 - 1 1 0 

V 36 ’ 72’ 


L = 


IT IT IT 
2’2^’6 


IT TT 


(A4) 


v=||/,o,o 


tt tt 




Z = 


(o,o,|). 


For the "H-l lattice we use the ORCC Brillouin zone 
from Ref. 94. 

171 <^, 171 (A5) 

These generalize quite simply to other even and odd-n 
harmonic honeycomb lattices. 

For the 2D lattice we take a rectangular Brillouin zone. 


27r /— 

0 < — - 171/73 - \ky\. 


(A 6 ) 


2. Constraining the ground state flux 

Because the mirror symmetries on the odd-n lattices 
do not pass through any lattice points we can use Lieb’s 
theorem^^ to constrain the ground state flux through the 
plaquettes through which the mirror plane passes. On 
the "H-l lattice this forces the hexagons as well as the 
symmetric 14-site plaquettes to carry 0-flux. The “tt-Aux 
state” fulfills these requirements.®® That flux configura¬ 
tion can be realized by changing the sign of on ev¬ 

ery other bridge bond or by threading tt flux through ev¬ 
ery rhombus in a projection along the c-direction. How¬ 
ever, the gauge illustrated in Fig. 1 (with a sign change 
on every other x and y bond in the oi and 02 directions 
on one ladder type) is most convenient because it accom¬ 
modates the TT flux state with the smallest possible unit 
cell (i.e. double the unit cell of the "H-l lattice). 


3. Brillouin zones 


The even-n lattices are face-centered cubics (Fddd), 
while the odd-n are base-centered (Cccm).®^’®®’®®’®® The 
2D honeycomb lattice has point group D^d, which is dis¬ 
tinct from the others. The three lattices considered here 
are representatives of these three space groups. 

For the H-0 lattice our choice of BZ is ORCFi from 
table 9 of Ref. 94. It can be parametrized by®^ 


|7| < 


TT 

3’ 


171 < 


77 


297r _ |7 _ 

36 72 3 ■ 

(A3) 


4. Hopping matrices 

The Majorana spinon Hamiltonian on the 'H-n lattice 
in the flux background $ is specified by the matrix 'D'^ 
from Eq. (10) (recall rj = n(4))). For each lattice the 
sites are numbered by their position in the positive c- 
direction. We choose a gauge such that = +1 if t 

is at an odd site (yellow) and j is at an even site (white) 


see Fig. 

1 . The hopping 

; matrices are®®’®® 




.j-,oo( 0 ) _ 

J" + 


g*ni.k 


< 

(A7) 



7 -, 0 ( 0 ) 

■^k 

= 

■ 

B- 

.r 


(AS) 



■ 

0 

0 





■^k “ 

AX 

"^k 

0 

K 

0 

0 

0 


(A9) 



0 

0 

B- 

J" 



II 

^ Ja + 

iki 


Sk 

= J| + 


K = 

: Jl + J3e- 

— iki 



= J% + 

1 


where = k • and the subscript on the couplings dis¬ 
tinguishes J® and jy bonds that lie on chains consisting 
of A and B bonds - see Eq. (9). The 7r-flux Hamiltonian 
requires a doubled unit cell. For a basis we take 

C (ci, C3, C5, C7, , ..., , C2, C4, Cg, Cg, C2, . •., Cg ) , 

(AlO) 
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where we have dropped the k-dependence to simplify the 
expression. In the gauge illustrated in Fig. 1 (the sign 
of is switched on the circled bonds), with the unit 
cell doubled in the 02 direction, the hopping is described 

by 

■^k ~ 

■ J" 0 0 0 0 0 0 ' 

.r 0 0 0 0 0 0 

0 Jl 0 0 0 0 

0 0 J% J" 0 0 0 

0 0 0 0 J" 0 0 ’ 

0 0 0 0 0 0 
0 0 0 0 Jl r 0 

0 0 0 0 0 J" _ 

where we define 

A^- =-J% + Al- = -Jl + . 

(All) 

We have confirmed that the energy of the tt-Aux state 
is lower on the "H-l lattice by studying the translation 
invariant system in the thermodynamic limit with an en¬ 
larged unit cell.®® In particular, ed, = J 

the 0-flux state has energy —0.77397J per site while the 
TT-flux has —0.77611J. For comparison, the local flux gap 
on the "H-O lattice is ^ O.IJ.^® 


Appendix B: Explicit Raman correlation fnnctions 

Here we show the equivalence between Eqs. (25) 
and (20) in the main text, provided that the Raman 
operator is of the form given in Eq. (18). We consider 
the response of a quadratic fermionic Hamiltonian in the 
presence of a quadratic fermionic Raman operator. As¬ 
suming translation invariance, Eqs. (25) and (18) can be 
written: 

I{uj) = J , 

R = ^ i?k, 

k 

Inserting the Fourier decomposition for R into the ex¬ 
pression for the Raman intensity gives 

/(a;)= /dt^e*‘"‘(0|i?k(OAk'(0)|0). (HI) 

k,k' 

Notice that only the Bjnm',k terms survive in (Bl) at 
zero temperature since the Amm',k terms annihilate the 


ground state. Therefore, 

/(cc) = ^ ( —S(uj £m,k ^n,k) 

k 

X (0| k®m',-kOm.k) (^^nn' ,k0^n,k^n' ,-kj I®) ) 

where we have used that ^(f) = jj.(0) in the 

Heisenberg picture. Applying the anti-commutation re¬ 
lations and writing the sum over bands explicitly we find: 

— TT ^ ^ ( (5 (ca £^m,k ,k^Rrnm' 

k^O m,m' 

(B2) 

The non-trivial commutation relations between oj^. and 
alk for k = —k lead a cancellation of the k = 0 term. 
Since the density of state vanishes at k = 0 the canceled 
term also vanishes in thermodynamic limit. This leads 
to Eq. (20) of the main text. 

Although the Majorana spinon Hamiltonian is 
translation-invariant in the flux sectors considered in the 
main text, the introduction of isolated fluxes would break 
this symmetry. To evaluate the Raman spectra at fi¬ 
nite flux density requires numerical simulations on finite¬ 
sized lattices. Though we have not carried these out, 
based on the numerics of Ref. 88 we expect that these 
translation-breaking terms will affect the qualitative con¬ 
clusions given in the main text only at low frequencies. 

Finally, a few numerical details. The Brillouin zone 
integral implicit in (20) is computed separately for each 
polarization combination. Numerically, we sample the 
Brilloiun zone and sort the points by energy. Then for a 
given w, the delta-function is replaced by a boxcar func¬ 
tion of finite width. Therefore, the integral is taken as a 
sum over points whose energies satisfy the delta function 
constraint up to a small tolerance in energy. A similar 
integration method was used in Ref. 38. Error estimates 
are obtained by randomly shifting the mesh and recom¬ 
puting the results. Typical error estimates for the value 
of I{uj) at a given value of oj are < 0.2%. 


Appendix C: Polarization dependence 

As discussed in the main text, the Raman operator 
R = X]a(^in ■ d“)(ein • d“)iJ“ is generally a different op¬ 
erator for different polarization choices of incoming and 
outgoing light given by ein and Cout, leading to many dif¬ 
ferent spectra denoted Here we discuss the 

relations between Raman intensities that come from sym¬ 
metry and the explicit form of the operators. 

If we act on the system with a space group symme¬ 
try transformation, the Raman intensity should be the 
same. That is, the symmetries must act trivially on the 
Raman intensity. We first consider the action of the sym¬ 
metries on the Raman operator. This corresponds to 
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changing the polarization so that for a symmetry action 
parametrized by we get 

R —> E V'(^out 

, 1 ^' 

In particular R^i, —)■ O^^iR^i^iO^'u. The action of the 
symmetries on the symmetric tensor R^^, must form a 
representation of the space group, which is a quadratic 
representation since it has two spatial indices. The task is 
then to associate the operators in the basis with the 
irreps. Then, since the symmetries must act trivially on 
/, the cross terms between R^^, associated with different 
irreps must vanish. 

1. 3D lattices 

For all n < oo the point group is D2h- Since this 
is orthorhombic there are no symmetries mixing the a, b, 
and c directions so that none of the with = a,b^c 
are mixed by the symmetry operations. This allows us 
to identify each of the 6 independent operators with an 
irrep. The irreps of D2h can be found in standard tables. 
It has eight ID irreps. However, since inversion can only 
act trivially on quadratic operators, R cannot contain 
any component in four of them. The four irreps that 
transform trivially under inversion are Aig, Big, B2g, and 
Bsg. 

Using Eq. (Cl) one can check that the aa, and cc 
components of R^^i, transform trivially under every point 
group symmetry operation and therefore belong to the 
trivial Aig irrep. The other three components, a6, ac, 
and be, transform under Big,B2g, and B^g respectively. 
In particular, two of the mirror symmetries act non- 
trivially in each of these irreps (e.g. Rat °> —Rat)- 
To summarize, the symmetries imply that the only non¬ 
zero intensity cross terms are Iaa,bb, Iaa,cc, and Ibb,cc and 
the other ones are zero (Recall Eq. (25)). 

Up to now we have only used the point group and 
thus these results hold for any Hamiltonian respecting 
the point group symmetries of the lattice. For a spin- 
exchange Hamiltonian we find further simplification upon 
writing out the Loudon-Fleury operator (Recall Eqs. (15) 
and (16)). For the bond vectors (9) we get 

^aa = H%+H% +Hi + Hi (C2) 

2V2Rab = HI - HI + Hi-Hi (C3) 

2V2Rb, = -HI + H%+Hl-Hl (C4) 

4Rac = -HI - HI + Hi + Hi (C5) 

Rcc = Raa + H\ (C6) 

First we find Rbb = 2Raa- This follows from the effective 
screw rotation discussed in the main text. Second, note 
that ZRaa + Rcc = H. We call this relationship, which 
comes from the Loudon-Fleury form of the Raman oper¬ 
ator, the LF relationship. 


Consider what this means for the correlation functions 
I. In the ground state = 0. Therefore, 

since Rcc |0) = —3Raa |0) + E |0), where E is the energy 
of |0) we find that 

{RgAt)RcM) = -3 (B^.(t)Baa(O)). (C7) 

This gives Ice = 9/aa = —3/aa,cc for any nearest-neighbor 
spin-exchange Hamiltonian on these lattices. 

The combination of symmetries and relations imposed 
by the explicit form of the operators Rgi, reduce the num¬ 
ber of independent, non-zero spectra to four. We choose 
to plot laai lab, lac, and I^c as representatives. 

2. 2D Honeycomb 

The point group of the honeycomb lattice is 
which has 3 irreps: Aig, A 2 g, and Eg. The Aig and 
A 2 g irreps are both ID but Eg is 2D because it involves 
nontrivial representation of Cg rotation. In A 2 g each of 
the three reflection is represented by —1, which cannot 
be realized by a quadratic operator. We say that this 
channel is not Raman active. 

Unlike the case for the 3D lattices, the Cg rotation 
symmetries do mix Raman operators in the Rgu basis 
with fj,,!/ = x,y. In particular, the linear combination 
RAig = transforms under the single irrep Aig 

and the vector {Rxy, Rsg) transforms under the Eg irrep, 
where Re = Qjie can check that the action of 

the symmetry elements on this vector realizes the full 2D 
representation of the Cg rotations given in Eq. (A2). 

Following the same procedure as for the 3D lattices we 
write out the Hamiltonian into its separate parts on each 
different bond (by orientation): H = J2a=x y z Then 
for any nearest-neighbor spin-exchange Hamiltonian re¬ 
specting the lattice symmetry the LF Raman operators 
have the form 

Rxx = 1{H^+ Ry), Ryy =H^ + ^Rxx, (C8) 

RA^g = H, Rsg = -H^ E^Rxx, (C9) 

Rxy="fiHy-Hn. (CIO) 

We see that the RAig channel does not create any exci¬ 
tations since it is equal to the Hamiltonian. Therefore 
only the Eg irrep is active on this lattice. 

We can relate Rxy and REg by considering the action of 
Cg rotation on the correlation functions involving those 
operators. We find 

1 3 V^3 

Ixy,xy —t -^Ixy,xy + -^lEg,Eg ± —^lEg,xy (CH) 

This implies that lEg,xy = 0 and that Ixy,xy = lEg,Eg- 

The fact that the Hamiltonian makes no excitations 
also implies that the operators H^ + H^ and —H^ give 
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rise to the same correlation functions and, by some rear¬ 
rangement, Rxx, —Ryy, and REg lead to the same Raman 
intensities. Finally, we have Ixx,xx = Iyy,yy = -Ixx,yy = 
Ixy,xy Combining these relationships with the symmetry 
constraints we have Ixx,xx — ^yy,yy ~ ^xy,xy — ^xx,yy 
and Ixx,xy = 0 = lyy^xy, which leaves only one inde¬ 
pendent spectrum. Considering an arbitrary polariza¬ 
tion Cin = cos 6 x + sin 6 y and Cout = cos (j)x. + sin (j)y we 
can compute I in terms of Ixx,xx, for instance. In this 
case one finds that the angular dependence drops out and 
I = Ixx,xx independent of polarization. 

When we take in the Kitaev Hamilto¬ 

nian we break the down to D 2 h- In fact, this group 
remains the point group even if we let ^ J'^. In this 
group the quadratic operators Rxx and Ryy act sepa¬ 
rately under the Aig irrep (they are even under reflec¬ 
tion) and Rxy falls under Big. The i? 2 g and B^g irreps 
are not active because there are no bonds along the di¬ 


rection out of the honeycomb plane. The LF relationship 
still gives Ixx,xx = lyy.yy = -Ixx,yy and the symmetry 
gives Ixy,xx = 0 = Ixy.yy, reducing I to the sum of two 
independent correlation functions. We choose to plot Ixx 
and Ixy . 

The inclusion of next-nearest-neighbors (NNNs) is 
quite simple under this formalism - one simply treats 
these as different bonds over which a runs. For the hon¬ 
eycomb lattice there is no expected change in the po¬ 
larization dependence because Rxx + Ryy = H holds 
for any bond orientation. However, on the 3D lat¬ 
tices we find that the relation to the Hamimltonian for 
NNNs cannot be satisfied simultaneously with the NN 
one 3i?aa + Rcc = H because the coefhcient 3 in front of 
Raa results from the specific of the bond-orientations that 
occur at nearest-neighbor only. Therefore the differences 
in Ice and 9Iaa, for instance, are a probe of terms beyond 
a NN spin-exchange Hamiltonian. 
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